# Load Necessary Packages
library(APCI)
library(ggplot2)
library(tidyr)
library(dplyr)
library(ggpubr)

# Load Datasets
setwd("/users/slajp8/Downloads/Monarchism in the UK")
data <- read.csv("IE data RR.csv",header=T,sep=",")

# Figure 1
cbbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
data.overtime.1 <- data %>% dplyr::group_by(period) %>% dplyr::summarise(m=weighted.mean(monarchism,weight,na.rm=T),a=weighted.mean(abolition,weight,na.rm=T),nhs=weighted.mean(nhsapproval,weight,na.rm=T),gov=weighted.mean(trustgov,weight,na.rm=T),mps=weighted.mean(trustmps,weight,na.rm=T))
data.overtime <- data.frame(period=rep(data.overtime.1$period,5),value=c(data.overtime.1$m,data.overtime.1$a,data.overtime.1$nhs,data.overtime.1$gov,data.overtime.1$mps),series=factor(c(rep("Monarchism",23),rep("Abolition",23),rep("NHS",23),rep("Govt",23),rep("MPs",23)),levels=c("Monarchism","Abolition","NHS","Govt","MPs")))
## FILL IN MISSING VALUES SOLELY FOR PLOTTING PURPOSES
data.overtime.filled <- data.frame(data.overtime)
data.overtime.filled$value[57] <- 0.49899538
data.overtime.filled$value[58] <- 0.49899538
data.overtime.filled$value[60] <- 0.56859741
data.overtime.filled$value[61] <- 0.56859741
data.overtime.filled$value[72] <- 0.35559899
data.overtime.filled$value[76] <- 0.37737122
data.overtime.filled$value[83] <- 0.36812993
data.overtime.filled$value[86] <- 0.29084007
data.overtime.filled$value[88] <- 0.30444107
data.overtime.filled$value[89] <- 0.30444107
data.overtime.filled$value[90] <- 0.30444107
data.overtime.filled$value[93] <- 0.22292615
data.overtime.filled$value[95] <- 0.19812655
data.overtime.filled$value[99] <- 0.19937565
data.overtime.filled$value[106] <- 0.19994452
data.overtime.filled$value[109] <- 0.17440927
data.overtime.filled$value[110] <- 0.17440927
data.overtime.filled$value[111] <- 0.17440927
data.overtime.filled$value[112] <- 0.17440927
data.overtime.filled$value[113] <- 0.17440927
p.overtime <- ggplot(data.overtime.filled,aes(x=period,y=value,group=series,color=series)) + geom_point() + geom_line() + theme_classic() + scale_color_manual(name="Variable",values=cbbPalette) + ylim(0,1) + ylab("Support") + xlab("Year")
# ggsave("Figure 1.png",p.overtime,device="png",width=6,height=4,units="in")

# Run Baseline APCs
data.80 <- subset(data,agegroup<=80)
data.80$agegroup <- factor(data.80$agegroup,levels=c("17","20","25","30","35","40","45","50","55","60","65","70","75","80"))
data.80$periodgroup <- factor(data.80$periodgroup,levels=c("1982","1992","1997","2002","2007","2012","2017","2022"))
## Monarchism
mod.baseline.m <- APCI::temp_model(data=data.80,outcome='monarchism',age='agegroup',period='periodgroup',covariate=c('country'),family='gaussian',weight='weight')
cohort.baseline.m <- APCI::cohortdeviation(mod.baseline.m$A,mod.baseline.m$P,mod.baseline.m$C,model=mod.baseline.m$model,covariate=c('country'),weight='weight')
## Abolition
mod.baseline.a <- APCI::temp_model(data=data.80,outcome='abolition',age='agegroup',period='periodgroup',covariate=c('country'),family='gaussian',weight='weight')
cohort.baseline.a <- APCI::cohortdeviation(mod.baseline.a$A,mod.baseline.a$P,mod.baseline.a$C,model=mod.baseline.a$model,covariate=c('country'),weight='weight')

# Figure 2
## Grand Means
weighted.mean(data$monarchism,w=data$weight,na.rm=T) # Monarchism=0.69
weighted.mean(data$abolition,w=data$weight,na.rm=T) # Abolition=0.09
## Monarchism
f2data.m <- data.frame(age=c(15,20,25,30,35,40,45,50,55,60,65,70,75,80),
                       effect=summary(mod.baseline.m$model)$coefficients[4:17,1],
                       se=summary(mod.baseline.m$model)$coefficients[4:17,2])
f2data.m$lci <- f2data.m$effect-(1.96*f2data.m$se)
f2data.m$uci <- f2data.m$effect+(1.96*f2data.m$se)
f2data.m$sig <- factor(ifelse(f2data.m$lci<0 & f2data.m$uci>0,"No","Yes"))
## Abolition
f2data.a <- data.frame(age=c(15,20,25,30,35,40,45,50,55,60,65,70,75,80),
                       effect=summary(mod.baseline.a$model)$coefficients[4:17,1],
                       se=summary(mod.baseline.a$model)$coefficients[4:17,2])
f2data.a$lci <- f2data.a$effect-(1.96*f2data.a$se)
f2data.a$uci <- f2data.a$effect+(1.96*f2data.a$se)
f2data.a$sig <- factor(ifelse(f2data.a$lci<0 & f2data.a$uci>0,"No","Yes"))
## Joint Figure
f2.m <- ggplot(f2data.m,aes(x=age,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Age") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Age, Monarchism") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.15, 0.17)
f2.a <- ggplot(f2data.a,aes(x=age,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Age") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Age, Abolition") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.07,0.09)
f2 <- ggpubr::ggarrange(f2.m,f2.a,nrow=1)
# ggsave("Figure 2.png",f2,device="png",width=6,height=4,units="in")

# Figure 3
## Monarchism
f3data.m <- data.frame(cohort=seq(1900,2000,by=5),
                       effect=as.numeric(as.character(cohort.baseline.m$cohort_average$cohort_average)),
                       lci=as.numeric(as.character(cohort.baseline.m$cohort_average$cohort_average))-(1.96*as.numeric(as.character(cohort.baseline.m$cohort_average$cohort_average_se))),
                       uci=as.numeric(as.character(cohort.baseline.m$cohort_average$cohort_average))+(1.96*as.numeric(as.character(cohort.baseline.m$cohort_average$cohort_average_se))))
f3data.m$sig <- ifelse(f3data.m$lci<0 & f3data.m$uci>0,"No","Yes")
## Abolition
f3data.a <- data.frame(cohort=seq(1900,2000,by=5),
                       effect=as.numeric(as.character(cohort.baseline.a$cohort_average$cohort_average)),
                       lci=as.numeric(as.character(cohort.baseline.a$cohort_average$cohort_average))-(1.96*as.numeric(as.character(cohort.baseline.a$cohort_average$cohort_average_se))),
                       uci=as.numeric(as.character(cohort.baseline.a$cohort_average$cohort_average))+(1.96*as.numeric(as.character(cohort.baseline.a$cohort_average$cohort_average_se))))
f3data.a$sig <- ifelse(f3data.a$lci<0 & f3data.a$uci>0,"No","Yes")
## Joint Figure
f3.m <- ggplot(f3data.m,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Cohort Average, Monarchism") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.20,0.16)
f3.a <- ggplot(f3data.a,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Cohort Average, Abolition") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.18, 0.13)
f3 <- ggpubr::ggarrange(f3.m,f3.a,nrow=1)
# ggsave("Figure 3.png",f3,device="png",width=6,height=4,units="in")

# Figure 4
## Monarchism
f4data.m <- data.frame(cohort=seq(1900,2000,by=5),
                       effect=as.numeric(as.character(cohort.baseline.m$cohort_slope$cohort_slope)),
                       lci=as.numeric(as.character(cohort.baseline.m$cohort_slope$cohort_slope))-(1.96*as.numeric(as.character(cohort.baseline.m$cohort_slope$cohort_slope_se))),
                       uci=as.numeric(as.character(cohort.baseline.m$cohort_slope$cohort_slope))+(1.96*as.numeric(as.character(cohort.baseline.m$cohort_slope$cohort_slope_se))))
f4data.m$sig <- ifelse(f4data.m$lci<0 & f4data.m$uci>0,"No","Yes")
## Abolition
f4data.a <- data.frame(cohort=seq(1900,2000,by=5),
                       effect=as.numeric(as.character(cohort.baseline.a$cohort_slope$cohort_slope)),
                       lci=as.numeric(as.character(cohort.baseline.a$cohort_slope$cohort_slope))-(1.96*as.numeric(as.character(cohort.baseline.a$cohort_slope$cohort_slope_se))),
                       uci=as.numeric(as.character(cohort.baseline.a$cohort_slope$cohort_slope))+(1.96*as.numeric(as.character(cohort.baseline.a$cohort_slope$cohort_slope_se))))
f4data.a$sig <- ifelse(f4data.a$lci<0 & f4data.a$uci>0,"No","Yes")
## Joint Figure
f4.m <- ggplot(f4data.m,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Cohort Slopes, Monarchism") + theme(plot.title = element_text(hjust = 0.5))
f4.a <- ggplot(f4data.a,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Cohort Slopes, Abolition") + theme(plot.title = element_text(hjust = 0.5))
f4 <- ggpubr::ggarrange(f4.m,f4.a,nrow=1)
# ggsave("Figure 4.png",f4,device="png",width=6,height=4,units="in")

# Figure A1
## Monarchism
fa1data.m <- data.frame(period=c(1992,1997,2002,2007,2012,2017,2022),
                        effect=summary(mod.baseline.m$model)$coefficients[19:25,1],
                        se=summary(mod.baseline.m$model)$coefficients[19:25,2])
fa1data.m$lci <- fa1data.m$effect-(1.96*fa1data.m$se)
fa1data.m$uci <- fa1data.m$effect+(1.96*fa1data.m$se)
fa1data.m$sig <- factor(ifelse(fa1data.m$lci<0 & fa1data.m$uci>0,"No","Yes"))
## Abolition
fa1data.a <- data.frame(period=c(1992,1997,2002,2007,2012,2017,2022),
                        effect=summary(mod.baseline.a$model)$coefficients[19:25,1],
                        se=summary(mod.baseline.a$model)$coefficients[19:25,2])
fa1data.a$lci <- fa1data.a$effect-(1.96*fa1data.a$se)
fa1data.a$uci <- fa1data.a$effect+(1.96*fa1data.a$se)
fa1data.a$sig <- factor(ifelse(fa1data.a$lci<0 & fa1data.a$uci>0,"No","Yes"))
## Joint Figure
fa1.m <- ggplot(fa1data.m,aes(x=period,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Period") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Period, Monarchism") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.06,0.1)
fa1.a <- ggplot(fa1data.a,aes(x=period,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Period") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Period, Abolition") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.08,0.06)
fa1 <- ggpubr::ggarrange(fa1.m,fa1.a,nrow=1)
# ggsave("Figure A1.png",fa1,device="png",width=6,height=4,units="in")

# APC of NHS Satisfaction
## Models
mod.nhs.m <- APCI::temp_model(data=data.80,outcome='nhsapproval',age='agegroup',period='periodgroup',covariate=c('country'),family='gaussian',weight='weight')
cohort.nhs.m <- APCI::cohortdeviation(mod.nhs.m$A,mod.nhs.m$P,mod.nhs.m$C,model=mod.nhs.m$model,covariate=c('country'),weight='weight')
## Age Plot
nhsdata.a <- data.frame(age=c(20,25,30,35,40,45,50,55,60,65,70,75,80),
                       effect=summary(mod.nhs.m$model)$coefficients[4:16,1],
                       se=summary(mod.nhs.m$model)$coefficients[4:16,2])
nhsdata.a$lci <- nhsdata.a$effect-(1.96*nhsdata.a$se)
nhsdata.a$uci <- nhsdata.a$effect+(1.96*nhsdata.a$se)
nhsdata.a$sig <- factor(ifelse(nhsdata.a$lci<0 & nhsdata.a$uci>0,"No","Yes"))
p.nhs.a <- ggplot(nhsdata.a,aes(x=age,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Age") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Age") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.1, 0.23)
## Period Plot
nhsdata.p <- data.frame(period=c(1992,1997,2002,2007,2012,2017,2022),
                        effect=summary(mod.nhs.m$model)$coefficients[17:23,1],
                        se=summary(mod.nhs.m$model)$coefficients[17:23,2])
nhsdata.p$lci <- nhsdata.p$effect-(1.96*nhsdata.p$se)
nhsdata.p$uci <- nhsdata.p$effect+(1.96*nhsdata.p$se)
nhsdata.p$sig <- factor(ifelse(nhsdata.p$lci<0 & nhsdata.p$uci>0,"No","Yes"))
p.nhs.p <- ggplot(nhsdata.p,aes(x=period,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Period") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Period") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.1, 0.23)
## Cohort Average
nhsdata.ca <- data.frame(cohort=seq(1900,2000,by=5),
                         effect=as.numeric(as.character(cohort.nhs.m$cohort_average$cohort_average)),
                         lci=as.numeric(as.character(cohort.nhs.m$cohort_average$cohort_average))-(1.96*as.numeric(as.character(cohort.nhs.m$cohort_average$cohort_average_se))),
                         uci=as.numeric(as.character(cohort.nhs.m$cohort_average$cohort_average))+(1.96*as.numeric(as.character(cohort.nhs.m$cohort_average$cohort_average_se))))
nhsdata.ca$sig <- ifelse(nhsdata.ca$lci<0 & nhsdata.ca$uci>0,"No","Yes")
p.nhs.ca <- ggplot(nhsdata.ca,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Cohort Average") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.1, 0.23)
## Cohort Slope
nhsdata.cs <- data.frame(cohort=seq(1900,2000,by=5),
                         effect=as.numeric(as.character(cohort.nhs.m$cohort_slope$cohort_slope)),
                         lci=as.numeric(as.character(cohort.nhs.m$cohort_slope$cohort_slope))-(1.96*as.numeric(as.character(cohort.nhs.m$cohort_slope$cohort_slope_se))),
                         uci=as.numeric(as.character(cohort.nhs.m$cohort_slope$cohort_slope))+(1.96*as.numeric(as.character(cohort.nhs.m$cohort_slope$cohort_slope_se))))
nhsdata.cs$sig <- ifelse(nhsdata.cs$lci<0 & nhsdata.cs$uci>0,"No","Yes")
p.nhs.cs <- ggplot(nhsdata.cs,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Cohort Slope") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.1, 0.23)
## Together
p.nhs <- ggpubr::ggarrange(p.nhs.a,p.nhs.p,p.nhs.ca,p.nhs.cs,nrow=2,ncol=2)
# ggsave("NHS APCI.png",p.nhs,device="png",width=6,height=4,units="in")

# APC of Political Trust
## Models
data.trust <- subset(data.80,periodgroup!="1982")
data.trust$agegroup[data.trust$agegroup=="17"] <- "20"
data.trust$periodgroup <- factor(data.trust$periodgroup,levels=c("1992","1997","2002","2007","2012","2017","2022"))
mod.govtrust.m <- APCI::temp_model(data=data.trust,outcome='trustgov',age='agegroup',period='periodgroup',covariate=c('country'),family='gaussian',weight='weight')
cohort.govtrust.m <- APCI::cohortdeviation(mod.govtrust.m$A,mod.govtrust.m$P,mod.govtrust.m$C,model=mod.govtrust.m$model,covariate=c('country'),weight='weight')
## Age Plot
trustdata.a <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80),
                        effect=summary(mod.govtrust.m$model)$coefficients[4:15,1],
                        se=summary(mod.govtrust.m$model)$coefficients[4:15,2])
trustdata.a$lci <- trustdata.a$effect-(1.96*trustdata.a$se)
trustdata.a$uci <- trustdata.a$effect+(1.96*trustdata.a$se)
trustdata.a$sig <- factor(ifelse(trustdata.a$lci<0 & trustdata.a$uci>0,"No","Yes"))
p.trust.a <- ggplot(trustdata.a,aes(x=age,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Age") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Age") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.15, 0.09)
## Period Plot
trustdata.p <- data.frame(period=c(1997,2002,2007,2012,2017,2022),
                        effect=summary(mod.govtrust.m$model)$coefficients[16:21,1],
                        se=summary(mod.govtrust.m$model)$coefficients[16:21,2])
trustdata.p$lci <- trustdata.p$effect-(1.96*trustdata.p$se)
trustdata.p$uci <- trustdata.p$effect+(1.96*trustdata.p$se)
trustdata.p$sig <- factor(ifelse(trustdata.p$lci<0 & trustdata.p$uci>0,"No","Yes"))
p.trust.p <- ggplot(trustdata.p,aes(x=period,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Period") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Period") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.15, 0.09)
## Cohort Average
trustdata.ca <- data.frame(cohort=seq(1910,2000,by=5),
                         effect=as.numeric(as.character(cohort.govtrust.m$cohort_average$cohort_average)),
                         lci=as.numeric(as.character(cohort.govtrust.m$cohort_average$cohort_average))-(1.96*as.numeric(as.character(cohort.govtrust.m$cohort_average$cohort_average_se))),
                         uci=as.numeric(as.character(cohort.govtrust.m$cohort_average$cohort_average))+(1.96*as.numeric(as.character(cohort.govtrust.m$cohort_average$cohort_average_se))))
trustdata.ca$sig <- ifelse(trustdata.ca$lci<0 & trustdata.ca$uci>0,"No","Yes")
p.trust.ca <- ggplot(trustdata.ca,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Cohort Average") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.15, 0.09)
## Cohort Slope
trustdata.cs <- data.frame(cohort=seq(1910,2000,by=5),
                         effect=as.numeric(as.character(cohort.govtrust.m$cohort_slope$cohort_slope)),
                         lci=as.numeric(as.character(cohort.govtrust.m$cohort_slope$cohort_slope))-(1.96*as.numeric(as.character(cohort.govtrust.m$cohort_slope$cohort_slope_se))),
                         uci=as.numeric(as.character(cohort.govtrust.m$cohort_slope$cohort_slope))+(1.96*as.numeric(as.character(cohort.govtrust.m$cohort_slope$cohort_slope_se))))
trustdata.cs$sig <- ifelse(trustdata.cs$lci<0 & trustdata.cs$uci>0,"No","Yes")
p.trust.cs <- ggplot(trustdata.cs,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Cohort Slope") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.15, 0.09)
## Together
p.trust <- ggpubr::ggarrange(p.trust.a,p.trust.p,p.trust.ca,p.trust.cs,nrow=2,ncol=2)
# ggsave("Gov Trust APCI.png",p.trust,device="png",width=6,height=4,units="in")

# APC of MP Trust
## Models
data.mptrust <- subset(data.trust,periodgroup!="2017")
mod.mpstrust.m <- APCI::temp_model(data=data.mptrust,outcome='trustmps',age='agegroup',period='periodgroup',covariate=c('country'),family='gaussian',weight='weight')
cohort.mpstrust.m <- APCI::cohortdeviation(mod.mpstrust.m$A,mod.mpstrust.m$P,mod.mpstrust.m$C,model=mod.mpstrust.m$model,covariate=c('country'),weight='weight')
## Age Plot
mpdata.a <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80),
                          effect=summary(mod.mpstrust.m$model)$coefficients[4:15,1],
                          se=summary(mod.mpstrust.m$model)$coefficients[4:15,2])
mpdata.a$lci <- mpdata.a$effect-(1.96*mpdata.a$se)
mpdata.a$uci <- mpdata.a$effect+(1.96*mpdata.a$se)
mpdata.a$sig <- factor(ifelse(mpdata.a$lci<0 & mpdata.a$uci>0,"No","Yes"))
p.mps.a <- ggplot(mpdata.a,aes(x=age,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Age") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Age") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.1, 0.09)
## Period Plot
mpdata.p <- data.frame(period=c(2002,2007,2012,2017,2022),
                          effect=summary(mod.mpstrust.m$model)$coefficients[16:20,1],
                          se=summary(mod.mpstrust.m$model)$coefficients[16:20,2])
mpdata.p$lci <- mpdata.p$effect-(1.96*mpdata.p$se)
mpdata.p$uci <- mpdata.p$effect+(1.96*mpdata.p$se)
mpdata.p$sig <- factor(ifelse(mpdata.p$lci<0 & mpdata.p$uci>0,"No","Yes"))
p.mps.p <- ggplot(mpdata.p,aes(x=period,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Period") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Period") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.1, 0.09)
## Cohort Average
mpdata.ca <- data.frame(cohort=seq(1915,2000,by=5),
                           effect=as.numeric(as.character(cohort.mpstrust.m$cohort_average$cohort_average)),
                           lci=as.numeric(as.character(cohort.mpstrust.m$cohort_average$cohort_average))-(1.96*as.numeric(as.character(cohort.mpstrust.m$cohort_average$cohort_average_se))),
                           uci=as.numeric(as.character(cohort.mpstrust.m$cohort_average$cohort_average))+(1.96*as.numeric(as.character(cohort.mpstrust.m$cohort_average$cohort_average_se))))
mpdata.ca$sig <- ifelse(mpdata.ca$lci<0 & mpdata.ca$uci>0,"No","Yes")
p.mps.ca <- ggplot(mpdata.ca,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Cohort Average") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.1, 0.09)
## Cohort Slope
mpdata.cs <- data.frame(cohort=seq(1915,2000,by=5),
                           effect=as.numeric(as.character(cohort.mpstrust.m$cohort_slope$cohort_slope)),
                           lci=as.numeric(as.character(cohort.mpstrust.m$cohort_slope$cohort_slope))-(1.96*as.numeric(as.character(cohort.mpstrust.m$cohort_slope$cohort_slope_se))),
                           uci=as.numeric(as.character(cohort.mpstrust.m$cohort_slope$cohort_slope))+(1.96*as.numeric(as.character(cohort.mpstrust.m$cohort_slope$cohort_slope_se))))
mpdata.cs$sig <- ifelse(mpdata.cs$lci<0 & mpdata.cs$uci>0,"No","Yes")
p.mps.cs <- ggplot(mpdata.cs,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Cohort Slope") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.1, 0.09)
## Together
p.mps <- ggpubr::ggarrange(p.mps.a,p.mps.p,p.mps.ca,p.mps.cs,nrow=2,ncol=2)
# ggsave("MP Trust APCI.png",p.mps,device="png",width=6,height=4,units="in")

# APC controlling for socio-demographics
## Models
mod.controls.m <- APCI::temp_model(data=data.trust,outcome='monarchism',age='agegroup',period='periodgroup',covariate=c('country','educ','class','female'),family='gaussian',weight='weight')
cohort.controls.m <- APCI::cohortdeviation(mod.controls.m$A,mod.controls.m$P,mod.controls.m$C,model=mod.controls.m$model,covariate=c('country','educ','class','female'),weight='weight')
mod.controls.a <- APCI::temp_model(data=data.trust,outcome='abolition',age='agegroup',period='periodgroup',covariate=c('country','educ','class','female'),family='gaussian',weight='weight')
cohort.controls.a <- APCI::cohortdeviation(mod.controls.a$A,mod.controls.a$P,mod.controls.a$C,model=mod.controls.a$model,covariate=c('country','educ','class','female'),weight='weight')
## Age Plots
controldata.m.a <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80),
                       effect=summary(mod.controls.m$model)$coefficients[7:18,1],
                       se=summary(mod.controls.m$model)$coefficients[7:18,2])
controldata.m.a$lci <- controldata.m.a$effect-(1.96*controldata.m.a$se)
controldata.m.a$uci <- controldata.m.a$effect+(1.96*controldata.m.a$se)
controldata.m.a$sig <- factor(ifelse(controldata.m.a$lci<0 & controldata.m.a$uci>0,"No","Yes"))
controldata.a.a <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80),
                              effect=summary(mod.controls.a$model)$coefficients[7:18,1],
                              se=summary(mod.controls.a$model)$coefficients[7:18,2])
controldata.a.a$lci <- controldata.a.a$effect-(1.96*controldata.a.a$se)
controldata.a.a$uci <- controldata.a.a$effect+(1.96*controldata.a.a$se)
controldata.a.a$sig <- factor(ifelse(controldata.a.a$lci<0 & controldata.a.a$uci>0,"No","Yes"))
p.control.m.a <- ggplot(controldata.m.a,aes(x=age,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Age") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Black"),guide="none") + ggtitle("Demogs") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.15, 0.17)
p.control.a.a <- ggplot(controldata.a.a,aes(x=age,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Age") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Demogs") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.07,0.09)
## Period Plots
controldata.m.p <- data.frame(period=c(1997,2002,2007,2012,2017,2022),
                          effect=summary(mod.controls.m$model)$coefficients[19:24,1],
                          se=summary(mod.controls.m$model)$coefficients[19:24,2])
controldata.m.p$lci <- controldata.m.p$effect-(1.96*controldata.m.p$se)
controldata.m.p$uci <- controldata.m.p$effect+(1.96*controldata.m.p$se)
controldata.m.p$sig <- factor(ifelse(controldata.m.p$lci<0 & controldata.m.p$uci>0,"No","Yes"))
controldata.a.p <- data.frame(period=c(1997,2002,2007,2012,2017,2022),
                              effect=summary(mod.controls.a$model)$coefficients[19:24,1],
                              se=summary(mod.controls.a$model)$coefficients[19:24,2])
controldata.a.p$lci <- controldata.a.p$effect-(1.96*controldata.a.p$se)
controldata.a.p$uci <- controldata.a.p$effect+(1.96*controldata.a.p$se)
controldata.a.p$sig <- factor(ifelse(controldata.a.p$lci<0 & controldata.a.p$uci>0,"No","Yes"))
p.control.m.p <- ggplot(controldata.m.p,aes(x=period,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Period") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Demogs") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.06,0.1)
p.control.a.p <- ggplot(controldata.a.p,aes(x=period,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Period") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Demogs") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.08,0.06)
## Cohort Average
controldata.m.ca <- data.frame(cohort=seq(1910,2000,by=5),
                           effect=as.numeric(as.character(cohort.controls.m$cohort_average$cohort_average)),
                           lci=as.numeric(as.character(cohort.controls.m$cohort_average$cohort_average))-(1.96*as.numeric(as.character(cohort.controls.m$cohort_average$cohort_average_se))),
                           uci=as.numeric(as.character(cohort.controls.m$cohort_average$cohort_average))+(1.96*as.numeric(as.character(cohort.controls.m$cohort_average$cohort_average_se))))
controldata.m.ca$sig <- ifelse(controldata.m.ca$lci<0 & controldata.m.ca$uci>0,"No","Yes")
controldata.a.ca <- data.frame(cohort=seq(1910,2000,by=5),
                               effect=as.numeric(as.character(cohort.controls.a$cohort_average$cohort_average)),
                               lci=as.numeric(as.character(cohort.controls.a$cohort_average$cohort_average))-(1.96*as.numeric(as.character(cohort.controls.a$cohort_average$cohort_average_se))),
                               uci=as.numeric(as.character(cohort.controls.a$cohort_average$cohort_average))+(1.96*as.numeric(as.character(cohort.controls.a$cohort_average$cohort_average_se))))
controldata.a.ca$sig <- ifelse(controldata.a.ca$lci<0 & controldata.a.ca$uci>0,"No","Yes")
p.control.m.ca <- ggplot(controldata.m.ca,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Demogs") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.20,0.16)
p.control.a.ca <- ggplot(controldata.a.ca,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Demogs") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.18, 0.13)
## Cohort Slope
controldata.m.cs <- data.frame(cohort=seq(1910,2000,by=5),
                           effect=as.numeric(as.character(cohort.controls.m$cohort_slope$cohort_slope)),
                           lci=as.numeric(as.character(cohort.controls.m$cohort_slope$cohort_slope))-(1.96*as.numeric(as.character(cohort.controls.m$cohort_slope$cohort_slope_se))),
                           uci=as.numeric(as.character(cohort.controls.m$cohort_slope$cohort_slope))+(1.96*as.numeric(as.character(cohort.controls.m$cohort_slope$cohort_slope_se))))
controldata.m.cs$sig <- ifelse(controldata.m.cs$lci<0 & controldata.m.cs$uci>0,"No","Yes")
controldata.a.cs <- data.frame(cohort=seq(1910,2000,by=5),
                               effect=as.numeric(as.character(cohort.controls.a$cohort_slope$cohort_slope)),
                               lci=as.numeric(as.character(cohort.controls.a$cohort_slope$cohort_slope))-(1.96*as.numeric(as.character(cohort.controls.a$cohort_slope$cohort_slope_se))),
                               uci=as.numeric(as.character(cohort.controls.a$cohort_slope$cohort_slope))+(1.96*as.numeric(as.character(cohort.controls.a$cohort_slope$cohort_slope_se))))
controldata.a.cs$sig <- ifelse(controldata.a.cs$lci<0 & controldata.a.cs$uci>0,"No","Yes")
p.control.m.cs <- ggplot(controldata.m.cs,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Demogs") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.13, 0.11)
p.control.a.cs <- ggplot(controldata.a.cs,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Demogs") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.11, 0.17)

# APC controlling for socio-demographics + ideology (LR, AL)
## Models
mod.controls1.m <- APCI::temp_model(data=data.trust,outcome='monarchism',age='agegroup',period='periodgroup',covariate=c('country','educ','class','female','lr','al'),family='gaussian',weight='weight')
cohort.controls1.m <- APCI::cohortdeviation(mod.controls1.m$A,mod.controls1.m$P,mod.controls1.m$C,model=mod.controls1.m$model,covariate=c('country','educ','class','female','lr','al'),weight='weight')
mod.controls1.a <- APCI::temp_model(data=data.trust,outcome='abolition',age='agegroup',period='periodgroup',covariate=c('country','educ','class','female','lr','al'),family='gaussian',weight='weight')
cohort.controls1.a <- APCI::cohortdeviation(mod.controls1.a$A,mod.controls1.a$P,mod.controls1.a$C,model=mod.controls1.a$model,covariate=c('country','educ','class','female','lr','al'),weight='weight')
## Age Plots
controldata1.m.a <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80),
                              effect=summary(mod.controls1.m$model)$coefficients[9:20,1],
                              se=summary(mod.controls1.m$model)$coefficients[9:20,2])
controldata1.m.a$lci <- controldata1.m.a$effect-(1.96*controldata1.m.a$se)
controldata1.m.a$uci <- controldata1.m.a$effect+(1.96*controldata1.m.a$se)
controldata1.m.a$sig <- factor(ifelse(controldata1.m.a$lci<0 & controldata1.m.a$uci>0,"No","Yes"))
controldata1.a.a <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80),
                              effect=summary(mod.controls1.a$model)$coefficients[9:20,1],
                              se=summary(mod.controls1.a$model)$coefficients[9:20,2])
controldata1.a.a$lci <- controldata1.a.a$effect-(1.96*controldata1.a.a$se)
controldata1.a.a$uci <- controldata1.a.a$effect+(1.96*controldata1.a.a$se)
controldata1.a.a$sig <- factor(ifelse(controldata1.a.a$lci<0 & controldata1.a.a$uci>0,"No","Yes"))
p.control1.m.a <- ggplot(controldata1.m.a,aes(x=age,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Age") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("D + Ideology") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.15, 0.17)
p.control1.a.a <- ggplot(controldata1.a.a,aes(x=age,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Age") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("D + Ideology") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.07,0.09)
## Period Plots
controldata1.m.p <- data.frame(period=c(1997,2002,2007,2012,2017,2022),
                              effect=summary(mod.controls1.m$model)$coefficients[21:26,1],
                              se=summary(mod.controls1.m$model)$coefficients[21:26,2])
controldata1.m.p$lci <- controldata1.m.p$effect-(1.96*controldata1.m.p$se)
controldata1.m.p$uci <- controldata1.m.p$effect+(1.96*controldata1.m.p$se)
controldata1.m.p$sig <- factor(ifelse(controldata1.m.p$lci<0 & controldata1.m.p$uci>0,"No","Yes"))
controldata1.a.p <- data.frame(period=c(1997,2002,2007,2012,2017,2022),
                              effect=summary(mod.controls1.a$model)$coefficients[21:26,1],
                              se=summary(mod.controls1.a$model)$coefficients[21:26,2])
controldata1.a.p$lci <- controldata1.a.p$effect-(1.96*controldata1.a.p$se)
controldata1.a.p$uci <- controldata1.a.p$effect+(1.96*controldata1.a.p$se)
controldata1.a.p$sig <- factor(ifelse(controldata1.a.p$lci<0 & controldata1.a.p$uci>0,"No","Yes"))
p.control1.m.p <- ggplot(controldata1.m.p,aes(x=period,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Period") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("D + Ideology") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.06,0.1)
p.control1.a.p <- ggplot(controldata1.a.p,aes(x=period,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Period") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("D + Ideology") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.08,0.06)
## Cohort Average
controldata1.m.ca <- data.frame(cohort=seq(1910,2000,by=5),
                             effect=as.numeric(as.character(cohort.controls1.m$cohort_average$cohort_average)),
                             lci=as.numeric(as.character(cohort.controls1.m$cohort_average$cohort_average))-(1.96*as.numeric(as.character(cohort.controls1.m$cohort_average$cohort_average_se))),
                             uci=as.numeric(as.character(cohort.controls1.m$cohort_average$cohort_average))+(1.96*as.numeric(as.character(cohort.controls1.m$cohort_average$cohort_average_se))))
controldata1.m.ca$sig <- ifelse(controldata1.m.ca$lci<0 & controldata1.m.ca$uci>0,"No","Yes")
controldata1.a.ca <- data.frame(cohort=seq(1910,2000,by=5),
                                effect=as.numeric(as.character(cohort.controls1.a$cohort_average$cohort_average)),
                                lci=as.numeric(as.character(cohort.controls1.a$cohort_average$cohort_average))-(1.96*as.numeric(as.character(cohort.controls1.a$cohort_average$cohort_average_se))),
                                uci=as.numeric(as.character(cohort.controls1.a$cohort_average$cohort_average))+(1.96*as.numeric(as.character(cohort.controls1.a$cohort_average$cohort_average_se))))
controldata1.a.ca$sig <- ifelse(controldata1.a.ca$lci<0 & controldata1.a.ca$uci>0,"No","Yes")
p.control1.m.ca <- ggplot(controldata1.m.ca,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("D + Ideology") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.20,0.16)
p.control1.a.ca <- ggplot(controldata1.a.ca,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("D + Ideology") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.18, 0.13)
## Cohort Slope
controldata1.m.cs <- data.frame(cohort=seq(1910,2000,by=5),
                             effect=as.numeric(as.character(cohort.controls1.m$cohort_slope$cohort_slope)),
                             lci=as.numeric(as.character(cohort.controls1.m$cohort_slope$cohort_slope))-(1.96*as.numeric(as.character(cohort.controls1.m$cohort_slope$cohort_slope_se))),
                             uci=as.numeric(as.character(cohort.controls1.m$cohort_slope$cohort_slope))+(1.96*as.numeric(as.character(cohort.controls1.m$cohort_slope$cohort_slope_se))))
controldata1.m.cs$sig <- ifelse(controldata1.m.cs$lci<0 & controldata1.m.cs$uci>0,"No","Yes")
controldata1.a.cs <- data.frame(cohort=seq(1910,2000,by=5),
                                effect=as.numeric(as.character(cohort.controls1.a$cohort_slope$cohort_slope)),
                                lci=as.numeric(as.character(cohort.controls1.a$cohort_slope$cohort_slope))-(1.96*as.numeric(as.character(cohort.controls1.a$cohort_slope$cohort_slope_se))),
                                uci=as.numeric(as.character(cohort.controls1.a$cohort_slope$cohort_slope))+(1.96*as.numeric(as.character(cohort.controls1.a$cohort_slope$cohort_slope_se))))
controldata1.a.cs$sig <- ifelse(controldata1.a.cs$lci<0 & controldata1.a.cs$uci>0,"No","Yes")
p.control1.m.cs <- ggplot(controldata1.m.cs,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("D + Ideology") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.13, 0.11)
p.control1.a.cs <- ggplot(controldata1.a.cs,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("D + Ideology") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.11, 0.17)

# APC controlling for socio-demographics + institutional trust
## Models
mod.controls2.m <- APCI::temp_model(data=data.trust,outcome='monarchism',age='agegroup',period='periodgroup',covariate=c('country','educ','class','female','trustgov'),family='gaussian',weight='weight')
cohort.controls2.m <- APCI::cohortdeviation(mod.controls2.m$A,mod.controls2.m$P,mod.controls3.m$C,model=mod.controls2.m$model,covariate=c('country','educ','class','female','trustgov'),weight='weight')
mod.controls2.a <- APCI::temp_model(data=data.trust,outcome='abolition',age='agegroup',period='periodgroup',covariate=c('country','educ','class','female','trustgov'),family='gaussian',weight='weight')
cohort.controls2.a <- APCI::cohortdeviation(mod.controls2.a$A,mod.controls2.a$P,mod.controls3.a$C,model=mod.controls2.a$model,covariate=c('country','educ','class','female','trustgov'),weight='weight')
## Age Plots
controldata2.m.a <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80),
                               effect=summary(mod.controls2.m$model)$coefficients[8:19,1],
                               se=summary(mod.controls2.m$model)$coefficients[8:19,2])
controldata2.m.a$lci <- controldata2.m.a$effect-(1.96*controldata2.m.a$se)
controldata2.m.a$uci <- controldata2.m.a$effect+(1.96*controldata2.m.a$se)
controldata2.m.a$sig <- factor(ifelse(controldata2.m.a$lci<0 & controldata2.m.a$uci>0,"No","Yes"))
controldata2.a.a <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80),
                               effect=summary(mod.controls2.a$model)$coefficients[8:19,1],
                               se=summary(mod.controls2.a$model)$coefficients[8:19,2])
controldata2.a.a$lci <- controldata2.a.a$effect-(1.96*controldata2.a.a$se)
controldata2.a.a$uci <- controldata2.a.a$effect+(1.96*controldata2.a.a$se)
controldata2.a.a$sig <- factor(ifelse(controldata2.a.a$lci<0 & controldata2.a.a$uci>0,"No","Yes"))
p.control2.m.a <- ggplot(controldata2.m.a,aes(x=age,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Age") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("D + Trust") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.15, 0.17)
p.control2.a.a <- ggplot(controldata2.a.a,aes(x=age,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Age") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("D + Trust") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.07,0.09)
## Period Plots
controldata2.m.p <- data.frame(period=c(1997,2002,2007,2012,2017,2022),
                               effect=summary(mod.controls2.m$model)$coefficients[20:25,1],
                               se=summary(mod.controls2.m$model)$coefficients[20:25,2])
controldata2.m.p$lci <- controldata2.m.p$effect-(1.96*controldata2.m.p$se)
controldata2.m.p$uci <- controldata2.m.p$effect+(1.96*controldata2.m.p$se)
controldata2.m.p$sig <- factor(ifelse(controldata2.m.p$lci<0 & controldata2.m.p$uci>0,"No","Yes"))
controldata2.a.p <- data.frame(period=c(1997,2002,2007,2012,2017,2022),
                               effect=summary(mod.controls2.a$model)$coefficients[20:25,1],
                               se=summary(mod.controls2.a$model)$coefficients[20:25,2])
controldata2.a.p$lci <- controldata2.a.p$effect-(1.96*controldata2.a.p$se)
controldata2.a.p$uci <- controldata2.a.p$effect+(1.96*controldata2.a.p$se)
controldata2.a.p$sig <- factor(ifelse(controldata2.a.p$lci<0 & controldata2.a.p$uci>0,"No","Yes"))
p.control2.m.p <- ggplot(controldata2.m.p,aes(x=period,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Period") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("D + Trust") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.06,0.1)
p.control2.a.p <- ggplot(controldata2.a.p,aes(x=period,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Period") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("D + Trust") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.08,0.06)
## Cohort Average
controldata2.m.ca <- data.frame(cohort=seq(1910,2000,by=5),
                              effect=as.numeric(as.character(cohort.controls2.m$cohort_average$cohort_average)),
                              lci=as.numeric(as.character(cohort.controls2.m$cohort_average$cohort_average))-(1.96*as.numeric(as.character(cohort.controls2.m$cohort_average$cohort_average_se))),
                              uci=as.numeric(as.character(cohort.controls2.m$cohort_average$cohort_average))+(1.96*as.numeric(as.character(cohort.controls2.m$cohort_average$cohort_average_se))))
controldata2.m.ca$sig <- ifelse(controldata2.m.ca$lci<0 & controldata2.m.ca$uci>0,"No","Yes")
controldata2.a.ca <- data.frame(cohort=seq(1910,2000,by=5),
                                effect=as.numeric(as.character(cohort.controls2.a$cohort_average$cohort_average)),
                                lci=as.numeric(as.character(cohort.controls2.a$cohort_average$cohort_average))-(1.96*as.numeric(as.character(cohort.controls2.a$cohort_average$cohort_average_se))),
                                uci=as.numeric(as.character(cohort.controls2.a$cohort_average$cohort_average))+(1.96*as.numeric(as.character(cohort.controls2.a$cohort_average$cohort_average_se))))
controldata2.a.ca$sig <- ifelse(controldata2.a.ca$lci<0 & controldata2.a.ca$uci>0,"No","Yes")
p.control2.m.ca <- ggplot(controldata2.m.ca,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("D + Trust") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.20,0.16)
p.control2.a.ca <- ggplot(controldata2.a.ca,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("D + Trust") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.18, 0.13)
## Cohort Slope
controldata2.m.cs <- data.frame(cohort=seq(1910,2000,by=5),
                              effect=as.numeric(as.character(cohort.controls2.m$cohort_slope$cohort_slope)),
                              lci=as.numeric(as.character(cohort.controls2.m$cohort_slope$cohort_slope))-(1.96*as.numeric(as.character(cohort.controls2.m$cohort_slope$cohort_slope_se))),
                              uci=as.numeric(as.character(cohort.controls2.m$cohort_slope$cohort_slope))+(1.96*as.numeric(as.character(cohort.controls2.m$cohort_slope$cohort_slope_se))))
controldata2.m.cs$sig <- ifelse(controldata2.m.cs$lci<0 & controldata2.m.cs$uci>0,"No","Yes")
controldata2.a.cs <- data.frame(cohort=seq(1910,2000,by=5),
                                effect=as.numeric(as.character(cohort.controls2.a$cohort_slope$cohort_slope)),
                                lci=as.numeric(as.character(cohort.controls2.a$cohort_slope$cohort_slope))-(1.96*as.numeric(as.character(cohort.controls2.a$cohort_slope$cohort_slope_se))),
                                uci=as.numeric(as.character(cohort.controls2.a$cohort_slope$cohort_slope))+(1.96*as.numeric(as.character(cohort.controls2.a$cohort_slope$cohort_slope_se))))
controldata2.a.cs$sig <- ifelse(controldata2.a.cs$lci<0 & controldata2.a.cs$uci>0,"No","Yes")
p.control2.m.cs <- ggplot(controldata2.m.cs,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("D + Trust") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.13, 0.11)
p.control2.a.cs <- ggplot(controldata2.a.cs,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("D + Trust") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.11, 0.17)

# APC controlling for socio-demographics + ideology + trust
## Models
mod.controls3.m <- APCI::temp_model(data=data.trust,outcome='monarchism',age='agegroup',period='periodgroup',covariate=c('country','educ','class','female','lr','al','trustgov'),family='gaussian',weight='weight')
cohort.controls3.m <- APCI::cohortdeviation(mod.controls3.m$A,mod.controls3.m$P,mod.controls3.m$C,model=mod.controls3.m$model,covariate=c('country','educ','class','female','lr','al','trustgov'),weight='weight')
mod.controls3.a <- APCI::temp_model(data=data.trust,outcome='abolition',age='agegroup',period='periodgroup',covariate=c('country','educ','class','female','lr','al','trustgov'),family='gaussian',weight='weight')
cohort.controls3.a <- APCI::cohortdeviation(mod.controls3.a$A,mod.controls3.a$P,mod.controls3.a$C,model=mod.controls3.a$model,covariate=c('country','educ','class','female','lr','al','trustgov'),weight='weight')
## Age Plots
controldata3.m.a <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80),
                               effect=summary(mod.controls3.m$model)$coefficients[10:21,1],
                               se=summary(mod.controls3.m$model)$coefficients[10:21,2])
controldata3.m.a$lci <- controldata3.m.a$effect-(1.96*controldata3.m.a$se)
controldata3.m.a$uci <- controldata3.m.a$effect+(1.96*controldata3.m.a$se)
controldata3.m.a$sig <- factor(ifelse(controldata3.m.a$lci<0 & controldata3.m.a$uci>0,"No","Yes"))
controldata3.a.a <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80),
                               effect=summary(mod.controls3.a$model)$coefficients[10:21,1],
                               se=summary(mod.controls3.a$model)$coefficients[10:21,2])
controldata3.a.a$lci <- controldata3.a.a$effect-(1.96*controldata3.a.a$se)
controldata3.a.a$uci <- controldata3.a.a$effect+(1.96*controldata3.a.a$se)
controldata3.a.a$sig <- factor(ifelse(controldata3.a.a$lci<0 & controldata3.a.a$uci>0,"No","Yes"))
p.control3.m.a <- ggplot(controldata3.m.a,aes(x=age,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Age") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("All") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.15, 0.17)
p.control3.a.a <- ggplot(controldata3.a.a,aes(x=age,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Age") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("All") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.07,0.09)
## Period Plots
controldata3.m.p <- data.frame(period=c(1997,2002,2007,2012,2017,2022),
                               effect=summary(mod.controls3.m$model)$coefficients[22:27,1],
                               se=summary(mod.controls3.m$model)$coefficients[22:27,2])
controldata3.m.p$lci <- controldata3.m.p$effect-(1.96*controldata3.m.p$se)
controldata3.m.p$uci <- controldata3.m.p$effect+(1.96*controldata3.m.p$se)
controldata3.m.p$sig <- factor(ifelse(controldata3.m.p$lci<0 & controldata3.m.p$uci>0,"No","Yes"))
controldata3.a.p <- data.frame(period=c(1997,2002,2007,2012,2017,2022),
                               effect=summary(mod.controls3.a$model)$coefficients[22:27,1],
                               se=summary(mod.controls3.a$model)$coefficients[22:27,2])
controldata3.a.p$lci <- controldata3.a.p$effect-(1.96*controldata3.a.p$se)
controldata3.a.p$uci <- controldata3.a.p$effect+(1.96*controldata3.a.p$se)
controldata3.a.p$sig <- factor(ifelse(controldata3.a.p$lci<0 & controldata3.a.p$uci>0,"No","Yes"))
p.control3.m.p <- ggplot(controldata3.m.p,aes(x=period,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Period") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("All") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.06,0.1)
p.control3.a.p <- ggplot(controldata3.a.p,aes(x=period,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Period") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("All") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.08,0.06)
## Cohort Average
controldata3.m.ca <- data.frame(cohort=seq(1910,2000,by=5),
                              effect=as.numeric(as.character(cohort.controls3.m$cohort_average$cohort_average)),
                              lci=as.numeric(as.character(cohort.controls3.m$cohort_average$cohort_average))-(1.96*as.numeric(as.character(cohort.controls3.m$cohort_average$cohort_average_se))),
                              uci=as.numeric(as.character(cohort.controls3.m$cohort_average$cohort_average))+(1.96*as.numeric(as.character(cohort.controls3.m$cohort_average$cohort_average_se))))
controldata3.m.ca$sig <- ifelse(controldata3.m.ca$lci<0 & controldata3.m.ca$uci>0,"No","Yes")
controldata3.a.ca <- data.frame(cohort=seq(1910,2000,by=5),
                                effect=as.numeric(as.character(cohort.controls3.a$cohort_average$cohort_average)),
                                lci=as.numeric(as.character(cohort.controls3.a$cohort_average$cohort_average))-(1.96*as.numeric(as.character(cohort.controls3.a$cohort_average$cohort_average_se))),
                                uci=as.numeric(as.character(cohort.controls3.a$cohort_average$cohort_average))+(1.96*as.numeric(as.character(cohort.controls3.a$cohort_average$cohort_average_se))))
controldata3.a.ca$sig <- ifelse(controldata3.a.ca$lci<0 & controldata3.a.ca$uci>0,"No","Yes")
p.control3.m.ca <- ggplot(controldata3.m.ca,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("All") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.20,0.16)
p.control3.a.ca <- ggplot(controldata3.a.ca,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("All") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.18, 0.13)
## Cohort Slope
controldata3.m.cs <- data.frame(cohort=seq(1910,2000,by=5),
                              effect=as.numeric(as.character(cohort.controls3.m$cohort_slope$cohort_slope)),
                              lci=as.numeric(as.character(cohort.controls3.m$cohort_slope$cohort_slope))-(1.96*as.numeric(as.character(cohort.controls3.m$cohort_slope$cohort_slope_se))),
                              uci=as.numeric(as.character(cohort.controls3.m$cohort_slope$cohort_slope))+(1.96*as.numeric(as.character(cohort.controls3.m$cohort_slope$cohort_slope_se))))
controldata3.m.cs$sig <- ifelse(controldata3.m.cs$lci<0 & controldata3.m.cs$uci>0,"No","Yes")
controldata3.a.cs <- data.frame(cohort=seq(1910,2000,by=5),
                                effect=as.numeric(as.character(cohort.controls3.a$cohort_slope$cohort_slope)),
                                lci=as.numeric(as.character(cohort.controls3.a$cohort_slope$cohort_slope))-(1.96*as.numeric(as.character(cohort.controls3.a$cohort_slope$cohort_slope_se))),
                                uci=as.numeric(as.character(cohort.controls3.a$cohort_slope$cohort_slope))+(1.96*as.numeric(as.character(cohort.controls3.a$cohort_slope$cohort_slope_se))))
controldata3.a.cs$sig <- ifelse(controldata3.a.cs$lci<0 & controldata3.a.cs$uci>0,"No","Yes")
p.control3.m.cs <- ggplot(controldata3.m.cs,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("All") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.13, 0.11)
p.control3.a.cs <- ggplot(controldata3.a.cs,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("All") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.11, 0.17)

# Figures Controlling
## Age
f2.m1 <- ggplot(f2data.m,aes(x=age,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Age") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("No Controls") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.15, 0.17)
f2.a1 <- ggplot(f2data.a,aes(x=age,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Age") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("No Controls") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.07,0.09)
p.m.age <- ggpubr::ggarrange(f2.m1,p.control.m.a,p.control1.m.a,p.control2.m.a,p.control3.m.a,nrow=2,ncol=3)
p.a.age <- ggpubr::ggarrange(f2.a1,p.control.a.a,p.control1.a.a,p.control2.a.a,p.control3.a.a,nrow=2,ncol=3)
# ggsave("Age Controlling Monarchism APCI.png",p.m.age,device="png",width=6,height=4,units="in")
# ggsave("Age Controlling Abolition APCI.png",p.a.age,device="png",width=6,height=4,units="in")
## Period
fa1.m1 <- ggplot(fa1data.m,aes(x=period,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Period") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("No Controls") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.06,0.1)
fa1.a1 <- ggplot(fa1data.a,aes(x=period,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Period") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("No Controls") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.08,0.06)
p.m.period <- ggpubr::ggarrange(fa1.m1,p.control.m.p,p.control1.m.p,p.control2.m.p,p.control3.m.p,nrow=2,ncol=3)
p.a.period <- ggpubr::ggarrange(fa1.a1,p.control.a.p,p.control1.a.p,p.control2.a.p,p.control3.a.p,nrow=2,ncol=3)
# ggsave("Period Controlling Monarchism APCI.png",p.m.period,device="png",width=6,height=4,units="in")
# ggsave("Period Controlling Abolition APCI.png",p.a.period,device="png",width=6,height=4,units="in")
## Cohort Averages
f3.m1 <- ggplot(f3data.m,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("No Controls") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.20,0.16)
f3.a1 <- ggplot(f3data.a,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("No Controls") + theme(plot.title = element_text(hjust = 0.5)) + ylim(-0.18, 0.13)
p.m.cohortavg <- ggpubr::ggarrange(f3.m1,p.control.m.ca,p.control1.m.ca,p.control2.m.ca,p.control3.m.ca,nrow=2,ncol=3)
p.a.cohortavg <- ggpubr::ggarrange(f3.a1,p.control.a.ca,p.control1.a.ca,p.control2.a.ca,p.control3.a.ca,nrow=2,ncol=3)
# ggsave("Cohort Average Controlling Monarchism APCI.png",p.m.cohortavg,device="png",width=6,height=4,units="in")
# ggsave("Cohort Average Controlling Abolition APCI.png",p.a.cohortavg,device="png",width=6,height=4,units="in")
## Cohort Slopes
f4.m1 <- ggplot(f4data.m,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("No Controls") + theme(plot.title = element_text(hjust = 0.5))
f4.a1 <- ggplot(f4data.a,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("No Controls") + theme(plot.title = element_text(hjust = 0.5))
p.m.cohortslope <- ggpubr::ggarrange(f4.m1,p.control.m.cs,p.control1.m.cs,p.control2.m.cs,p.control3.m.cs,nrow=2,ncol=3)
p.a.cohortslope <- ggpubr::ggarrange(f4.a1,p.control.a.cs,p.control1.a.cs,p.control2.a.cs,p.control3.a.cs,nrow=2,ncol=3)
# ggsave("Cohort Slope Controlling Monarchism APCI.png",p.m.cohortslope,device="png",width=6,height=4,units="in")
# ggsave("Cohort Slope Controlling Abolition APCI.png",p.a.cohortslope,device="png",width=6,height=4,units="in")

# Correlation between support for the monarchy and trust in other institutions
## NHS Approval
cor.test(subset(data,period==1983)$monarchism,subset(data,period==1983)$nhsapproval)
cor.test(subset(data,period==1994)$monarchism,subset(data,period==1994)$nhsapproval)
cor.test(subset(data,period==1995)$monarchism,subset(data,period==1995)$nhsapproval)
cor.test(subset(data,period==1996)$monarchism,subset(data,period==1996)$nhsapproval)
cor.test(subset(data,period==1997)$monarchism,subset(data,period==1997)$nhsapproval)
cor.test(subset(data,period==1998)$monarchism,subset(data,period==1998)$nhsapproval)
cor.test(subset(data,period==1999)$monarchism,subset(data,period==1999)$nhsapproval)
cor.test(subset(data,period==2000)$monarchism,subset(data,period==2000)$nhsapproval)
cor.test(subset(data,period==2002)$monarchism,subset(data,period==2002)$nhsapproval)
cor.test(subset(data,period==2003)$monarchism,subset(data,period==2003)$nhsapproval)
cor.test(subset(data,period==2005)$monarchism,subset(data,period==2005)$nhsapproval)
cor.test(subset(data,period==2006)$monarchism,subset(data,period==2006)$nhsapproval)
cor.test(subset(data,period==2007)$monarchism,subset(data,period==2007)$nhsapproval)
cor.test(subset(data,period==2008)$monarchism,subset(data,period==2008)$nhsapproval)
cor.test(subset(data,period==2011)$monarchism,subset(data,period==2011)$nhsapproval)
cor.test(subset(data,period==2012)$monarchism,subset(data,period==2012)$nhsapproval)
cor.test(subset(data,period==2015)$monarchism,subset(data,period==2015)$nhsapproval)
cor.test(subset(data,period==2017)$monarchism,subset(data,period==2017)$nhsapproval)
cor.test(subset(data,period==2018)$monarchism,subset(data,period==2018)$nhsapproval)
cor.test(subset(data,period==2021)$monarchism,subset(data,period==2021)$nhsapproval)
cor.test(subset(data,period==2022)$monarchism,subset(data,period==2022)$nhsapproval)
cor.test(subset(data,period==2023)$monarchism,subset(data,period==2023)$nhsapproval)
cor.test(subset(data,period==2024)$monarchism,subset(data,period==2024)$nhsapproval)
## Trust in Government
cor.test(subset(data,period==1983)$monarchism,subset(data,period==1983)$trustgov)
cor.test(subset(data,period==1994)$monarchism,subset(data,period==1994)$trustgov)
cor.test(subset(data,period==1995)$monarchism,subset(data,period==1995)$trustgov)
cor.test(subset(data,period==1996)$monarchism,subset(data,period==1996)$trustgov)
cor.test(subset(data,period==1997)$monarchism,subset(data,period==1997)$trustgov)
cor.test(subset(data,period==1998)$monarchism,subset(data,period==1998)$trustgov)
cor.test(subset(data,period==1999)$monarchism,subset(data,period==1999)$trustgov)
cor.test(subset(data,period==2000)$monarchism,subset(data,period==2000)$trustgov)
cor.test(subset(data,period==2002)$monarchism,subset(data,period==2002)$trustgov)
cor.test(subset(data,period==2003)$monarchism,subset(data,period==2003)$trustgov)
cor.test(subset(data,period==2005)$monarchism,subset(data,period==2005)$trustgov)
cor.test(subset(data,period==2006)$monarchism,subset(data,period==2006)$trustgov)
cor.test(subset(data,period==2007)$monarchism,subset(data,period==2007)$trustgov)
cor.test(subset(data,period==2008)$monarchism,subset(data,period==2008)$trustgov)
cor.test(subset(data,period==2011)$monarchism,subset(data,period==2011)$trustgov)
cor.test(subset(data,period==2012)$monarchism,subset(data,period==2012)$trustgov)
cor.test(subset(data,period==2015)$monarchism,subset(data,period==2015)$trustgov)
cor.test(subset(data,period==2017)$monarchism,subset(data,period==2017)$trustgov)
cor.test(subset(data,period==2018)$monarchism,subset(data,period==2018)$trustgov)
cor.test(subset(data,period==2021)$monarchism,subset(data,period==2021)$trustgov)
cor.test(subset(data,period==2022)$monarchism,subset(data,period==2022)$trustgov)
cor.test(subset(data,period==2023)$monarchism,subset(data,period==2023)$trustgov)
cor.test(subset(data,period==2024)$monarchism,subset(data,period==2024)$trustgov)
## Trust in MPs
cor.test(subset(data,period==1983)$monarchism,subset(data,period==1983)$trustmps)
cor.test(subset(data,period==1994)$monarchism,subset(data,period==1994)$trustmps)
cor.test(subset(data,period==1995)$monarchism,subset(data,period==1995)$trustmps)
cor.test(subset(data,period==1996)$monarchism,subset(data,period==1996)$trustmps)
cor.test(subset(data,period==1997)$monarchism,subset(data,period==1997)$trustmps)
cor.test(subset(data,period==1998)$monarchism,subset(data,period==1998)$trustmps)
cor.test(subset(data,period==1999)$monarchism,subset(data,period==1999)$trustmps)
cor.test(subset(data,period==2000)$monarchism,subset(data,period==2000)$trustmps)
cor.test(subset(data,period==2002)$monarchism,subset(data,period==2002)$trustmps)
cor.test(subset(data,period==2003)$monarchism,subset(data,period==2003)$trustmps)
cor.test(subset(data,period==2005)$monarchism,subset(data,period==2005)$trustmps)
cor.test(subset(data,period==2006)$monarchism,subset(data,period==2006)$trustmps)
cor.test(subset(data,period==2007)$monarchism,subset(data,period==2007)$trustmps)
cor.test(subset(data,period==2008)$monarchism,subset(data,period==2008)$trustmps)
cor.test(subset(data,period==2011)$monarchism,subset(data,period==2011)$trustmps)
cor.test(subset(data,period==2012)$monarchism,subset(data,period==2012)$trustmps)
cor.test(subset(data,period==2015)$monarchism,subset(data,period==2015)$trustmps)
cor.test(subset(data,period==2017)$monarchism,subset(data,period==2017)$trustmps)
cor.test(subset(data,period==2018)$monarchism,subset(data,period==2018)$trustmps)
cor.test(subset(data,period==2021)$monarchism,subset(data,period==2021)$trustmps)
cor.test(subset(data,period==2022)$monarchism,subset(data,period==2022)$trustmps)
cor.test(subset(data,period==2023)$monarchism,subset(data,period==2023)$trustmps)
cor.test(subset(data,period==2024)$monarchism,subset(data,period==2024)$trustmps)